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Abstract 



We consider the problem of finding the ground state of a model type-II su- 
perconductor on the two-dimensional surface of a sphere, penetrated by N 
vortices. Numerical work shows the ground states to consist of a triangular 
network of the vortices with twelve five-coordinated centres. Values of N are 
found with particularly low energy ground states, due to structures of high 
symmetry. The large limit is treated within elasticity theory to compare 
with the triangular vortex lattice that forms the ground state on an infinite 
flat plane. Together with numerical work this demonstrates that the thermo- 
dynamic limit A/^ — > oo of the spherical system remains different from the flat 
plane due to the presence of twelve disclination defects. 
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I. INTRODUCTION 



The problem of constructing an optimum lattice-like structure over a curved surface has 
become an area of interest in diverse contexts within condensed matter physics. Some exam- 
ples are in work on flexible tethered membranes, Fullerene molecules, the Thomson problem 
of electrons on a spherical surface, and in models of two-dimensional systems using a spher- 
ical geometry to study both the Quantum Hall Effect (QHE)and thin film superconductors. 

The well known Fullerene molecules demonstrate how a low energy structure can be 
formed by folding a hexagonal lattice of carbon atoms (as is found in graphite) onto a closed 
surface, as long as twelve five-membered rings (pentagons) are present- a simple consequence 
of Euler's theorem. These pentagons are essentially disclination defects in the hexagonal 
lattice. The first discovered Fullerene molecule was Cgo in which each of the sixty atoms 
holds an identical symmetry position within the structure of the molecule, so the atoms 
reside on the surface of a sphere.0 This is a special case and as the number of carbon atoms 
increases in these molecules the shape may distort from a sphere to reduce the strain from the 
ideal hexagonal lattice over large areas. The interaction between the energy cost of bending 
the surface of the molecule and the strain energy within the molecule due to the disclination 
defects is of central importance in the question of the stability of different structures. [2-4] 

The same considerations are important in the behaviour of membranes with internal 
crystalline order and the ability to buckle out of the two-dimensional plane. A disclination 
defect may lower its energy by buckling. If this reduces the total energy of the defect 
such that it diverges at most logarithmically with the system size, then this raises the 
possibility of a buckling transition as defects begin to proliferate at some finite temperature.!! 
Alternatively, if we consider a closed membrane, such as a vesicle made from surfactant 
bilayers, then the interaction of the internal orientational order with the physical curvature 
may alter the shape,i or even the topology^ of the membrane. 

Some important and well studied models involve these problems but with the curvature 
of the system fixed. One instance is Thomson's problem of trying to find the lowest energy 
configuration for N electrons that are constrained to lie on the surface of a sphere. Although 
the problem was first proposed in the rather dated context of classical models of the atom,! 
it has been extensively studied recently. This is partly because of the general relevance 
to any physical system of this geometry, but also for its interest as an unsolved problem, 
providing a testing ground for various numerical optimization methods. [9-13] 

A different model where a two-dimensional system is restricted to a spherical surface 
has been studied both in the context of the QHE and in numerical studies on thin film 
superconductors. In this model, a magnetic field perpendicular to the surface is imposed 
by placing a Dirac monopole at the centre of the sphere. The possible electronic states on 
the surface of the sphere split into Landau levels under the influence of the magnetic field, 
in analogy with the problem on a fiat plane. The reason the model was put on a sphere by 
Haldane was to allow "the construction of homogeneous states" with only a finite number 
of electrons.Q Recent Monte Carlo simulations of 2D superconductors have been performed 
with a similar model to Haldane's with a superconducting wavefunction in the lowest Landau 
level on the sphere. [20-22] This wavefunction contains N zeros that correspond to vortices 
in the supercurrent, where N depends on the quantized strength of the monopole. The 
reason for using this model is again to enable translational invariance which is not possible 
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in a finite system on a fiat plane. 

Tlie Monte Carlo simulations on a sphere have led us to consider in detail the ground 
states of this model. While the problem of finding the ground state of an infinite type II 
superconductor penetrated by vortices was long ago solved and found to be the triangular 
vortex lattice,lii such a lattice cannot form on a sphere without the presence of twelve vor- 
tices with only five nearest neighbours. These twelve vortices are the centres of disclination 
defects. By considering the strains from the perfect triangular lattice, caused by the discli- 
nation defects, within elasticity theory, we have approximated the ground state energy in the 
large N limit using similar methods used in theoretical work on membranes and FuUerene 
molecules. We have also found numerically the ground states with finite N, using symmetry 
considerations to reach large system sizes. We find values of that give particularly low 
energies and this is explained. By extrapolating our numerical results to large we find a 
finite energy cost per vortex on the sphere compared to the infinite plane ground state that 
is consistent with our results from elasticity theory. 



II. FORMULATION 



Our model thin film superconductor consists of a spherical shell of superconducting ma- 
terial, thickness d, radius R and a monopole at the centre of the sphere that produces an 
integer multiple of flux quanta through the spherical surface. We ignore spatial fluctuations 
in the magnetic flux density, B, at the surface; the effective penetration depth for super- 
currents in thin films becomes arbitrarily large as d is reduced. We choose a cylindrically 
symmetric gauge consistent with this field with A = (A^,, Ag, A^) = {0,0, BRta,n9/2). We 
measure lengths in the units Im = (^o/Svri?)^/^ which if there are N quanta of flux gives 
R = If we may describe the properties of the superconductor by a complex order 

parameter ip{6, 0), then the Ginzburg-Landau free energy Hamiltonian will be given by 



d^r 



4 + ^^*D2^ 

2m 



(2.1) 



where = D*. D and D = —ih'V — 2eA. We diagonalize the operator by expanding 
lb in a basis of eigenfunctions of which form degenerate Landau levels. The degenerate 
set with the lowest eigenvalue may be filled by the orthonormal functions^ 



^m(6',0) 



hr.e'"'^ sin™ 



(6/2) cos^-"'{9/2) 



(2.2) 



with m = 0,N and hm = [{N + l)!/47ri?^m!(A^ — m)!]^''^. This is the lowest Landau level 
(LLL) and over a large range of fields and temperatures it is a good approximation to 
restrict t/- to the LLL, ilj{e,<p) = Q E Vm^m{0 , (p) . We set Q = {^okBT / pdEf/^ . With this 
restriction we can write the Hamiltonian in terms of the basis coefficients, which for < 
is given by 



H{aT,{um]) = kBTa^T {{um}) 



N 



N 



m=0 



p,q,r,s=0 



(2.3) 
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where Wp+q^q^r is given inS, ar = dQ"^ («(^) + eBfi/m) /ksT is the reduced temperature 
variable and we have scaled the coefficients as Vm = Um\o:T\^^'^- The quartic term in Equa- 
tion ( p.3| ) can be rewritten to 



N 2N 
m=0 p=0 



(2.4) 



where = 2ttN Ef=o B'/\2N -p + l,p+ l)hqhp_qQ{p -q)Q(N + q- p)UqUp^q, B{x, y) = 
T{x)T{y)/T{x + y) is the beta function and 9(g) is the Heaviside step function. 

The vortices in this system correspond to the zeros in ip{6,(j)). The phase of the order 
parameter changes by 2tt when a path is followed that encircles any zero once. This becomes 
clear if we make the projection 

C = tan(^/2)e^'^. (2.5) 

This gives the form = cos^(^/2) E™=o = cos^(^/2)/jv(C)- Therefore /^(C) is a 

holomorphic function of ( with simple zeros in the complex-(^ plane. It can always be 
written in the alternative product form /Ar(C) = Cnili(C ~ d)- C is an overall complex 
amplitude and {Q} are the vortex positions in the projection of Equation ( ^.5] ). The function 
iIj{9, (p) is equally well described by the set {um} or the set {C, Q}. There is no simple relation 
between the two, although numerical routines may be used to find the positions of the zeros 
for a given set of basis coefficients. Despite this lack of a simple relation from the basis 
coefficients to coordinates on the sphere, there are still some spatial transformations one 
can make using the {um} formalism. For instance rotation about the z-axis by an angle 7 
may be performed by the transformation Um — ^ Um^^"^"'- Reflection in the x-z plane results 
in Um — ^ Rotation by n about the x-axis occurs under the change Um — > UN-m- 
We are interested in the ground states of JF ({«„}). We write the Hamiltonian as 

^=-A + ^A^ (2.6) 

where A = J^UmU^ and (3a = (|'^o|^)/(|V^o|^) = J2 I^pP/^^ is the Abrikosov factor. This 
is minimized by A = —N/Pa to give J^min = —NEq = —N/2Pa, so minimizing JF is 
equivalent to minimizing (3A{{um})- The correct A is then given by a scale factor on the 
basis coefficients that does not alter Pa- 

The ground state of the LLL vortex system on an infinite plane is well known to be 
the triangular lattice^ which has /3a = Paa — 1-1596. With periodic boundary conditions, 
this is also the ground state on the finite systems used in other simulations. However, 
a perfect triangular lattice cannot form on a spherical surface. The closest configuration 
the vortices can make to an ideal lattice must contain twelve "disclinations" , i.e. twelve 
vortices that only have five nearest neighbours. In Section ^ we give our results for directly 
minimizing J-" {{um}) using a simple numerical method, but first we describe in Section |T| 
calculations using elasticity theory to give the finite energy cost that the spherical system 
will have as ^ 00 due to the twelve disclinations. 
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III. ELASTICITY THEORY 



For a lattice in a two-dimensional plane, the elastic energy cost of deformations from the 
perfect ground state lattice is given in the harmonic approximation by 

Fei = lJ d?r {2iiul + AmL) , (3.1) 

where Uij{r) = |((9jtij(r) + djUiiv)) is the elastic strain matrix. The displacement 
u(r) = (uj;, Uy) represents the deformation of the lattice from the point r to the point r + u. 
The elastic constants and A are related to the shear and bulk moduli by Cghear = and 
Chuik = /i + A. For the LLL ground state the bulk modulus is infinite; the vortex system is 
incompressible so u^k = 0. The shear modulus is given byil 

/X = 0.48 X -/ioff^ V ^J^2 • (3.2) 

The GL parameter k is the ratio of magnetic and superconducting correlation lengths, which 
diverges when we neglect the magnetic screening of super currents. In our approximation we 
can write the shear modulus as /i = 0.48/cBTQ;y/47r/?^/^ = 0.0659i?o/^m- 

From the elastic energy in Equation Hooke's law may be derived by minimizing F^i 
to give diOij = for the stress tensor 

aij = 2fiUij + XukkSij. (3.3) 

The zero divergence condition allows the reformulation of the problem in terms of the Airy 
stress functioned aij = eikejidkdix- (This is analogous to the vector potential that ensures 
zero divergence of magnetic fields.) 

In the presence of topological defects, such as the disclinations we are considering, the 
displacement field u(r) is multi-valued. A disclination is defined by the change in bond angle 
closed loop is followed. Encircling a five-fold disclination in a triangular 
lattice will increase by 27r/6. This results in the noncommutativity of the derivatives of t? 
at the centre of the disclination. Writing the strain field in terms of the Airv stress function 
results in the biharmonic equation that contains all of 2D elasticity theory:B 

= sir). (3.4) 

The density of disclinations is s(r) = J2a ■^^^(r— Tq) where a labels each defect and Sa = 2tc/6 
for a five-fold disclination. In Equation (|3.4|) Kq is the 2D Young's modulus, which in the 
LLL is 

Ko = = 4^ = 0.264Eo/t. (3.5) 

Of course, our problem is not on a fiat plane but on a sphere, so we must take into 
account the bending of the system out of the plane. In the large N limit the surface will 
be flat locally compared to lattice spacings. Over a small region, we may approximate 
the sphere as a plane with some small perpendicular deflection, /(r). For our purposes we 
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neglect the bending energy which will tend to a constant — independent on the system size — 
when integrated over the whole sphere. However, we will need to write the strain matrix as 
\{diUj + djUi + difdjf). This alters the biharmonic equation by adding an extra term, 



u. 



«J 2 

det((9j(9j/) ~ the Gaussian curvature to Equation ( p.4|) : 

^V\ = s{v)-K{v). (3.6) 

For a sphere the curvature is constant, K = 1/B?. We write x &s the superposition of 
twelve contributions corresponding to each disclination, x(r) = Y^^a=iXaij)- The solution 
to Equation ( |3.6| ) is found in 0. The elastic energy from Equation ( |3.1| ) can be written in 
terms of the Airy stress function as in Equation ([A6|) . We have calculated the energy cost 
for these twelve disclinations for different configurations. For any configuration the total 
energy scales with the surface area of the sphere, so there is a finite energy cost per vortex 
in the limit N oo. The minimum energy is found with the disclinations at the corners of 
an icosahedron. In this case we find an energy cost per vortex of 

5Eei = OMTKoR^/N = O.OOSlEo. (3.7) 

In Section ^ we will compare 6Eei with the energy 6Eei obtained by direct minimization of 
^{{um}) for finite A^. 



IV. NUMERICAL RESULTS 

We have used a simple quasi- Newton algorithm to find configurations {«„} that minimize 
the Abrikosov ratio PaH'iJ'ti})] this is equivalent to finding the minimum energy of the system. 
Clearly there are some transformations of under which (3^ is invariant. The energy of 
the system remains unchanged after a global phase change in the order parameter, or after a 
rotation of the whole system about some axis. These freedoms can be fixed by restrictions on 
the coefficients.lll Our results up to A^ = 200 are shown in Figure |l|. The presence of "magic 
numbers" , for which the ground state has a lower energy than nearby values, is clearly seen 
at A^ = 12, 32, 72, 132 and 192. To some extent we can explain the magic numbers from 
the expected symmetry of the most stable ground states. 

By finding the zeros of the ground states, we can look at the vortex configurations. As 
might be expected they make up a triangular network, but with twelve five-coordinated 
centres. The magic number states display icosahedral symmetry with the disclinations cor- 
responding to the corners of the icosahedron, as in Figure ^ In fact the structures appear 
to be projections on to the sphere of icosadeltahedra, which are polyhedra with identi- 
cal equilateral triangular faces and icosahedral symmetry.il These may be constructed by 
considering the number of triangular lattice vectors between neighbouring five-fold centres, 
labelled by the indices {h,k). 
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FIG. 1. The minimum values of the Abrikosov ratio for different system sizes. Note the 
low values for N = 12, 32, 72, 132 and 192. The inset shows the values for large N plotted against 
A^~i/2 with a fit to extrapolate to the — > oo limit. The value for an infinite flat plane, /J^jO, is 
shown for comparison. 



(a) (b) (c) 

FIG. 2. The numerically found ground states for three magic number cases: (a) N = 32, (b) 
N = 72, (c) N = 132. Only one hemisphere is shown; the dots represent the zeros of the order 
parameter and the five-fold centres are represented by larger dots, which in each case form the 
corners of an icosahedron. 



A simple geometrical argument shows that this divides each face of the icosahedron into T 
triangular faces with T = + hk + k^ . so T may be equal to 1, 3, 4, 7, 9, 12, 13, 16, 19.... This 
gives a total of 20T faces and lOT + 2 vertices to the icosadeltahedra. The magic numbers 
we find satisfy these conditions, the only remaining question being why some possible values 
of N are not so low in energy, eg. = 42 or 122. 
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Notice that the structures for N = 72 and 132 possess a chirahty. In fact this will 
always be the case for h ^ k (as long as neither is zero). This may be a factor in deciding 
the lowest energy structures as the complex conjugation of the coefficients is equivalent to a 
reflection. If a ground state has no chirality then any reflection of the state will be equivalent 
to a rotation. This would require that we could write the ground state with all coefficients 
having the same phase. It is unlikely that this combination would be effective in minimizing 
the complex interactions in the quartic term of the Hamiltonian especially as the system 
size increases. Therefore the chirality allows a greater variety of phase differences between 
the coefficients. Such subtleties may explain why the ground states for some non-chiral 
icosahedral numbers are not particularly stable, as with = 42, 122. 

Remarkably, the same structures that are found when we minimize our vortex system 
for some magic numbers are seen in nature in the form of the shells of certain viruses.@ In 
particular the structure shown in Figure |](c) for N = 132 is also observed in double and 
single shelled simian rotaviruses, in a left handed configuration.i3 That these similarities 
exist in such different systems suggests that some general principle exists for the criterion 
of the most stable structure. It is possibly related to a mathematical problem that also 
generates these structures. This is the "covering" problem: How may N equal overlapping 
circles (without gaps) cover a sphere so that the diameter of the circles is minimized?c3 The 
centres of these circles correspond to the vortices in our system. The alternative problem of 
maximizing the diameter of the circles apparently gives different solutions. 

The details of the magic number states becomes less important as becomes large, 



the limit treated in Section |T| within elasticity theory. As A^ increases the numerical 



minimization becomes less trivial as there is an increasing overlap between the basis states 



resulting in more complex phase interference (see Ref. 29 where numerical minimization on 
a plane in cylindrical coordinates resulted in rather high values of (3a)- Another possibly 
related problem is the growth in the number of metastable states with energies only slightly 
above the ground state. The numerical work on Thomson's problem has found that the 
number of metastable states grows exponentially with the system size.0 We may use our 
knowledge of the symmetry of the most stable configurations to reduce the number of free 
variables by an order of magnitude. We assume the icosahedral properties of the ground 
state, and choose a five-coordinated vortex at 6 = 0. The five-fold symmetry about the z- 
axis is imposed by setting = for n 7^ (5m -|- 1) where m is any integer. The icosahedron 
also has five two-fold axes of symmetry at right angles to each five-fold axis. This two- 
fold rotational symmetry will arise if we set Un = ujsf-n- We have performed the same 
minimization routine using these constraints for A^ = 10m + 2 up to A^ = 652. 

Our results for large A^ are shown in the inset to Figure |^. The data fits well to the 
form f3A = A + BN'^/^ with A ^ 1.1624 = /5ao(1 + 0.0024) and B ~ 0.0648. As A^ ^ 00, 
(3a does not seem to converge to the infinite plane value (3ao (contrary to the conclusions of 



Ref. |T5| from the minimizations of small system sizes). This extrapolation implies a finite 



energy cost per vortex on the sphere in the large A^ limit of 

5E = 0.0024£'o. (4.1) 



The difference between this and the result of Equation ( |3.7| ) from elasticity theory may be 
explained by the inadequacy of the harmonic approximation in this calculation. As there 
are large strains associated with the disclination defects, non-linear effects will be important 
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in determining the total energy cost. In their calculation of disclination defects on flat 
membranes Nelson and Seunga found the same mismatch between elasticity theory and 
numerical results, with elasticity overestimating the energy cost by a similar proportion. 

Within the approximations of Section |T| the correction for finite N is not predicted. Nu- 
merically the energy cost per vortex falls as 1 /i? at large N which implies a total contribution 
that grows proportionally to the radius. 

It must be stressed that the results of this section required no great numerical effort. 
More sophisticated optimization methods (eg. simulated annealing! and its generalizationsEl) 
may give greater confidence in whether or not the absolute minima have been found. More 
extensive work would also give results for larger A^. However, our use of symmetry has 
allowed us to do a great deal with just a simple routine. 

V. CONCLUSIONS 

The original motivation of this work was for the use of the numerical ground states 
in Monte Carlo simulations.0 We also wanted to investigate the differences between these 
ground states and the ground states on a fiat plane, for which other groups have performed 
simulations obtaining different results. The work in this paper shows that the vortex ground 
states on the sphere do not approach the ground state of the infinite fiat plane as — oo. 
The presence of the twelve disclinations remains important however large the sphere. 

This work may also be of interest in wider contexts: first, in its relation to other optimiza- 
tion problems of points on a sphere. Our particular system allows a use of symmetry that 
may not be so straightforward with position variables. This enables us to find approximate 
ground states at large N with quite unsophisticated numerical techniques. Our elasticity 
calculation may be relevant to the large limit of Thomson's problem, using properties of 
the Wigner lattice on a 2D infinite plane. This limit has been considered before, and pro- 
jection of the Wigner lattice onto a spherical surface was used to estimate the extra energy 
on the spherclll However, no consideration was taken of the required disclination defects. 
This paper also provides a numerical test of the accuracy of elasticity theory for curved 
membranes and disclination defects where the approximation of small deviations from the 
ideal lattice breaks down. Finally, the fact that the structures we see in the magic number 
ground states are the same structures seen in such a different field as virology suggests that 
these fascinating shapes are the result of some general optimization criteria. 
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APPENDIX A: DISCLINATION ON A SPHERE 

In this appendix we derive the contributions from each of the twelve disclinations to the 
the Airy Stress function x(r) = I^iLi Xa(r) and describe how this leads to the elastic energy 
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cost of these disclinations. From Equation ( p.6|) with K = and s = | Z]q=i ^^(r ~ ^a) 
we have: 



Kq 6 



1 



(Al) 



where r — Tq, = {9' (a) , (j)' (a)) , and 9'{a) and 0'(q;) are the polar and azimuthal angles with 
respect to the axis through the disclination. Using the symmetry about this axis, this may 
be integrated to give 



V^Xa = ^(lnsin(^'(a)/2) + A), 



(A2) 



with A a constant. For Xa to be well defined the integral of V^x« over the whole sphere 
must be zero, which means A = In 2 — 1. Integrating again leads to 



dXo 



09' (a) 



24 



tan(0'(a)/2)(21nsin(0'(a)/2) + 3(ln2 - 1)) 



(A3) 



Xc 



In cos {9'{a)/2) [2 In sin {9'{a)/2) + 3(ln 2 - 1)] 



+ 



/ 

Jo 



X 



■dx}. 



(A4) 



The integral in Equation (^^) cannot be written as a finite number of elementary functions.!^ 
From Equation (|3.1|) the elastic energy of the disclinations can be written 



el 



1 

2 



-rri^^xf T7—^ik^ji9kdi (d^xdjx) 



(A5) 
(A6) 



with the 2D poisson ratio a = A/(2/i + A) equal to unity in the LLL approximation. The 
second term in Equation (|A6|) only gives contributions on boundaries, and so is zero on the 
sphere. Therefore from Equation ( |A^ ) we can find the energies of different configurations of 
the disclinations on the sphere. 
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